home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byte0987.arc / KAREX.ARC / KAREX1.BAS next >
Encoding:
BASIC Source File  |  1980-01-01  |  5.5 KB  |  251 lines
  1.  
  2. 100 '---------------------------------------
  3. 101 '
  4. 102 '  KAREX1.BAS is a Microsoft BASIC
  5.           Release 5 program
  6. 103 '  that solves EXAMPLE 1 of the article
  7. 104 '
  8. 105 '        KARMARKAR'S ALGORITHM
  9. 106 '
  10. 107 '    by Andrew M. Rockett
  11.         and John C. Stevenson
  12. 108 '
  13. 109 'This program was written
  14.       by Andrew M. Rockett
  15. 110 '
  16. 111 '---------------------------------------
  17. 200 '
  18. 202 ' N is number of unknowns and K is the
  19.           number of equations
  20. 204 '
  21. 206 N = 3 : K = 1
  22. 208 '
  23. 210 K1 = K + 1 : K2 = 2*K1
  24. 212 DIM A0(N), XOLD(N), XNEW(N), CC(N),
  25.     CP(N),
  26. A(K,N), B(K1,N), B1(K1,K2),
  27.     B2(N,K1), B3(N,N)
  28. 214 '
  29. 216 ' CC is for the objective function
  30. 218 ' B1, B2 and B3 are used for the
  31.       computation of CP
  32. 220 ' R and C are "row" and "column" indices
  33. 222 '
  34. 224 ' Initially, set XNew = A0, the center
  35.       of simplex
  36. 226 '
  37. 228 FOR C = 1 TO N: A0(C) = 1/N : XNEW(C) =
  38.     A0(C)
  39.     :NEXT C
  40. 230 '
  41. 232 ' T is the tolerance
  42. 234 '
  43. 236 T = .001
  44. 238 '
  45. 240 ' ALPHA is usually set equal to 1/4 ...
  46. 242 '
  47. 244 ALPHA = .25
  48. 246 '
  49. 248 ITERATION = 0
  50. 250 '
  51. 252 ' Data for constraint matrix A
  52. 254 '
  53. 256 DATA  2, -3, 1
  54. 258 '
  55. 260 FOR R = 1 TO K:FOR C = 1 TO N:READ
  56.     A(R,C):NEXT
  57.     C:NEXT R
  58. 262 '
  59. 264 ' Data for objective function CC
  60. 266 '
  61. 268 DATA 3,  3,  -1
  62. 270 '
  63. 272 FOR C = 1 TO N : READ CC(C) : NEXT C
  64. 274 '
  65. 276 ' Set initial Value to value at center
  66.       of simplex ...
  67. 278 '
  68. 280 V = 0 : FOR C=1 TO N:V = V + CC(C)*A0(C):NEXT 
  69.     C: VNEW = V
  70. 282 '
  71. 284 ' Now we can begin the MAIN ITERATION
  72.        process ...
  73. 286 '
  74. 300 WHILE VNEW/V > T
  75. 301 '
  76. 302 PRINT USING "####";ITERATION;:
  77.     FOR C=1 TO N:PRINT USING
  78.     "###.####";XNEW(C); :
  79.     NEXT C :
  80.  PRINT USING
  81.    "####.#######";VNEW/V
  82. 303 '
  83. 304 ITERATION = ITERATION + 1
  84. 305 '
  85. 306 ' Put Xnew into Xold
  86. 307 '
  87. 308 FOR C = 1 TO N : XOLD(C) = XNEW(C) : NEXT C
  88. 309 '
  89. 310 ' Construct the matrix B
  90. 311 '
  91. 312 FOR R=1 TO K:FOR C=1 TO
  92.     N:B(R,C)=A(R,C)*XOLD(C):
  93.     NEXT C:NEXT R
  94. 313 FOR C = 1 TO N : B(K1,C) = 1 : NEXT C
  95. 314 '
  96. 315 ' Zero matrices to be used in
  97.       computations...
  98. 316 '
  99. 317 FOR R=1 TO K1 : FOR C=1 TO K2 :
  100.     B1(R,C)=0 :
  101.     NEXT C : NEXT R
  102. 318 FOR R=1 TO N  : FOR C=1 TO K1 :
  103.     B2(R,C)=0 :
  104.     NEXT C : NEXT R
  105. 319 FOR R=1 TO N  : FOR C=1 TO N  :
  106.     B3(R,C)=0 :
  107.     NEXT C : NEXT R
  108. 320 FOR C=1 TO N  : CP(C) = 0 : NEXT C
  109. 321 '
  110. 322 ' Find BBT and put in B1
  111. 323 '
  112. 324 FOR R = 1 TO K1 :
  113.      FOR C = 1 TO K1 :
  114.       FOR I = 1 TO N:B1(R,C)=B1(R,C)+B(R,I)*B(C,I):
  115.        NEXT I:
  116.      NEXT C :
  117.     NEXT R
  118. 325 '
  119. 326 ' Adjoin an identity matrix to BBT
  120. 327 '
  121. 328 FOR I = 1 TO K1 : B1(I,I+K1)=1 : NEXT I
  122. 329 '
  123. 330 ' Row reduce BBT|I 
  124. 331 '
  125. 332 FOR R = 1 TO K1
  126. 333  IF B1(R,R) <> 0 THEN 338
  127. 334    I = R + 1
  128. 335     IF I > K1 THEN PRINT"Error! BBT is
  129.         SINGULAR!": GOTO 400
  130. 336     IF B1(I,R) = 0 THEN I = I+1 : GOTO
  131.         335
  132. 337     FOR C = 1 TO K2 : SWAP
  133.         B1(R,C),B1(I,C) : 
  134.         NEXT C
  135. 338  FOR I = R+1 TO K1 :
  136.        Z = B1(I,R)/B1(R,R):
  137.        FOR C=1 TO K2:B1(I,C)=B1(I,C)
  138.        -Z*B1(R,C):
  139.         NEXT C:
  140.      NEXT I
  141. 339 NEXT R
  142. 340 '
  143. 341 ' Now back substitute to finish it ...
  144. 342 '
  145. 343 FOR R = K1 TO 2 STEP -1 :
  146.      FOR I = R-1 TO 1 STEP -1 :
  147.        Z = B1(I,R)/B1(R,R) :
  148.        FOR C = R TO K2:B1(I,C)=B1(I,C)-Z*B1(R,C):
  149.         NEXT C:
  150.      NEXT I :
  151.     NEXT R
  152. 344 '
  153. 345 Remember to make diagonal entries 1's
  154. 346 '
  155. 347 FOR R=1 TO K1 :
  156.      Z = B1(R,R) :
  157.      FOR C = 1 TO K2 : B1(R,C) = B1(R,C)/Z :
  158.      NEXT C :
  159.     NEXT R
  160. 348 '
  161. 349 ' BBT Inverse is now in B1 in columns
  162.       K1+1 to K2
  163. 350 '
  164. 351 ' Now multiply BBT Inverse by BT and put
  165.       in B2
  166. 352 '
  167. 353 FOR R = 1 TO N :
  168.      FOR C = 1 TO K1 :
  169.       FOR J = 1 TO K1:B2(R,C)=
  170.        B2(R,C)+B(J,R)*B1
  171.        (J,C+K1):
  172.       NEXT J:
  173.      NEXT C :
  174.     NEXT R
  175. 354 '
  176. 355 ' Take THAT and multiply by B and put in
  177.       B3
  178. 356 '
  179. 357 FOR R = 1 TO N :
  180.      FOR C = 1 TO N :
  181.       FOR J = 1 TO
  182.        K1:B3(R,C)=B3(R,C)+B2(R,J)*B(J,C):
  183.        NEXT J:
  184.      NEXT C :
  185.     NEXT R
  186. 358 '
  187. 359 ' Find  I-B3 by subtracting 1's on
  188.        diagonal and changing signs
  189. 361 '
  190. 362 FOR R = 1 TO N : B3(R,R) = B3(R,R) - 1 :
  191.     NEXT R
  192. 363 FOR R=1 TO N:FOR C=1 TO N:B3(R,C) = 
  193.     -1*B3(R,C):
  194.     NEXT C:NEXT R
  195. 364 '
  196. 365 ' Multiply by D
  197. 366 '
  198. 367 FOR R=1 TO N:FOR C=1 TO
  199.     N:B3(R,C)=B3(R,C)*
  200.     XOLD(C):NEXT C:NEXT R
  201. 368 '
  202. 369 ' Find projection of CC and call it CP
  203. 370 '
  204. 371 FOR R=1 TO N:FOR C=1 TO
  205.     N:CP(R)=CP(R)+B3(R,C)*
  206.     CC(C):NEXT C:NEXT R
  207. 372 '
  208. 373 ' Find length of CP and the normalized
  209.       CP
  210. 374 '
  211. 375 AA = 0
  212. 376 FOR C=1 TO N : AA = AA + CP(C)*CP(C) :
  213.     NEXT C
  214. 377 AA = SQR(AA) : FOR C=1 TO N : CP(C) =
  215.     CP(C)/AA : NEXT C
  216. 378 '
  217. 379 'Find a*, project back to get new X ...
  218. 380 '
  219. 381 AA = SQR(N*(N-1))/ALPHA
  220. 382 FOR C=1 TO N : XNEW(C) = (A0(C) -
  221.     CP(C)/AA)*XOLD(C) : 
  222.     NEXT C
  223. 383 '
  224. 384 ' And remember to divide by "size" of
  225.       new X to complete the 385 ' projective
  226.  transformation back to the original simplex
  227. 386 '
  228. 387 AA = 0
  229. 388 FOR C=1 TO N : AA = AA + XNEW(C) : NEXT
  230.     C
  231. 389 FOR C=1 TO N : XNEW(C) = XNEW(C)/AA :
  232.     NEXT C
  233. 390 '
  234. 391 ' Find objective function Value at NEW
  235.       point X
  236. 392 '
  237. 393 VNEW = 0
  238. 394 FOR C=1 TO N : VNEW = VNEW + CC(C)*XNEW(C) : NEXT C
  239. 395 '
  240. 396 WEND '  End of main iteration loop ...
  241. 397 '
  242. 398 PRINT:PRINT"Tolerance reached:
  243.     Vnew/Vinitial = "; VNEW/V:PRINT
  244. 399 PRINT USING "####"; ITERATION; :
  245.     FOR C=1 TO N:PRINT USING
  246.         "###.####";XNEW(C); : 
  247.     NEXT C :
  248. PRINT USING
  249.        "####.#######";VNEW/V
  250. 400 END
  251.